#pragma once

#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <chrono>
#include <ctime>
#include <limits>
#include <vector>
#include <complex>
#include <cmath>
#include <Eigen/Dense>
#include <boost/math/special_functions/bessel.hpp>
#include <mkl.h>
#include <omp.h>

#define Complex std::complex<double>

enum Grid_type
{
	StaggeredTartanGrids = 1,
	TartanGrids = 2,
	ChebyshevGrids = 3,
	NewChebyshevGrids = 4
};

enum RBF_TYPE
{
	IMQ = 1,
	Bessel = 2,
	Cubic = 3,
	Quintic = 4,
	TPS = 5
};

class RBF
{
public:
	RBF(){};
	~RBF(){};

	std::pair<double, double> New_test_scattered_RBF(const double lengthx, const double lengthy, const double lengthz,
		const double lambda_, int gridtype,
		int Nums1Dx, int Nums1Dy, int Nums1Dz, const double c_new,
		std::vector<std::pair<std::vector<int>, std::vector<double>>>&
		Nums1Dxyz_Lengthxyz_OptimalC_list, int RBFtype);

private:
	std::vector<Eigen::Vector3d> r;
	Eigen::Vector3d rs;
	int Nums_1D_x, Nums_1D_y, Nums_1D_z;
	size_t nums_sp;
	double length_x, length_y, length_z, lambda, _k;
	std::vector<double> zeroX, zeroY, zeroZ, spx, spy, spz;
	std::vector<std::complex<double>> GG;
	std::vector<double> phi_mat;
	std::vector<std::vector<double>> TestXYZ_list;
	int RBF_type;
	int grid_type;

	void Init(const double lengthx, const double lengthy, const double lengthz, const double lambda_,
	          int gridtype, int Nums1Dx, int Nums1Dy, int Nums1Dz, int RBFtype);

	void generate_tartan_interpolation_zeros();

	void generate_Chebyshev_interpolation_zeros();

	void generate_NewChebyshev_interpolation_zeros();

	void generate_interpolation_zeros();

	double factorial(int n);

	void generate_GG();

	void generate_test_points();

	void generate_the_Monomials_matrix_fortran_format(const double c);

	void scattered_point_RBF_interpolation(Complex& Ga, double xx, double yy, double zz, const std::vector<std::vector<double>>& Phi_mat_Pinv, double c);

	double LagrangeBaseFunction(int idx, int k, double* tk, double t);

	void scattered_point_lagrange_interpolation(Complex& Ga, double xx, double yy, double zz);

	std::vector<std::vector<double>> Eigen_Cal_PseudoInverse();

	std::vector<std::vector<double>> MKL_Cal_PseudoInverse();

	void printCurrentTime(std::string printinfo);

};


